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Abstract 

Magnetic plasma instabilities appear to play an important role in the early stages of quark-gluon 
plasma equilibration in the high energy (weak coupling) limit. Numerical studies of the growth 
of such instabilities from small seed fluctuations have found initial exponential growth in their 
energy, followed by linear growth once the associated color magnetic fields become so large that 
their non-abelian interaction are non-perturbative. In this paper, we use simulations to determine 
the nature of this linear energy growth. We find that the long-wavelength modes associated with 
the instability have ceased to grow, but that they cascade energy towards the ultraviolet in the 
form of plasmon excitations of ever increasing energy. We find a quasi-steady-state power-law 
distribution oc k~ v for this cascade, with spectral index v ~ 2. 



I. INTRODUCTION 



One of the surprises which has emerged from the heavy ion experiments at RHIC is 
that the medium or plasma which is produced in a heavy ion collision appears to display 
collective behavior analogous to a fluid with a very small viscosity (the so-called "Quark- 
Gluon Liquid") [1]. In particular, hydrodynamic treatments taking the medium to be an 
ideal fluid give the correct flow description [2]. To clarify, we remind the reader that "ideal 
fluid" is the opposite of "ideal gas"; it means that scatterings or interactions so efficiently 
maintain local equilibrium (or at least isotropy), that the stress tensor remains everywhere 
isotropic when measured in the local rest frame. Gases behave like ideal fluids on distance 
and time scales large compared to transport mean free paths and mean free times. 

The most popular, and perhaps most likely, explanation for the small viscosity is that 
the Quark-Gluon Plasma at the energy densities achieved at RHIC is very far from weak 
coupling. The strong coupling a s is large, interactions are very strong, and the collective 
behavior is natural. To secure this interpretation, we would need to see that the observed 
behavior really differs from the expected weakly coupled behavior. Most treatments to date 
of weakly coupled hot QCD show that equilibration is slow [3] (see however [4]). However, 
we have recently shown [5] that even the most complete of these, the "bottom-up" scenario 
of Baier, Mueller, Schiff and Son [6], is incorrect, because it ignores plasma instabilities [7], 
which in fact dominate the physics of weakly-coupled anisotropic plasmas, in QCD or in 
ordinary QED. 

Plasma instabilities arise as a result of two pieces of physics. First, Lorentz contraction 
of the nuclei means that the initial region of plasma is a flat pancake shape. It is reasonable 
to expect at weak coupling that quarks and gluons will for a time fly in nearly straight lines 
at near the speed of light (at least, once the system has expanded enough that the densities 
of quarks and gluons are perturbative). Consider particles scattered in all directions in 
the initial moments of the collision. The starting geometry dictates that the momentum 
distribution of particles will subsequently become highly anisotropic with time, as illustrated 
in Fig. 1. Second, in the presence of such anisotropic momentum distributions, certain 
soft gauge field configurations grow exponentially, at least if they are small enough for a 
perturbative treatment to be reliable. 

Not only have plasma instabilities not been fully taken into account in studying the 
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FIG. 1: Cartoon of why the momentum distribution becomes anisotropic. The starting region is 
thin along the beam axis (vertical in this figure). Since particles fly in a straight line at the speed 
of light, particles at point * at time t originated at time on a sphere of radius t centered at *. 
But only part of this sphere was in the plasma; so only particles with small beam-axis momenta 
can get to point * at time t, and the momentum distribution is anisotropic. 
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evolution of the Quark-Gluon Plasma; instabilities in Yang-Mills theory in general are not 
yet well understood or well characterized. Therefore, at the current juncture we cannot 
say whether the experimental results at RHIC are in contradiction with weakly coupled 
predictions, because we simply do not know what the predictions of isotropization in a weakly 
coupled Quark-Gluon Plasma actually are. This rather embarrassing gap in our knowledge 
needs to be filled. In particular, the behavior of plasma instabilities in a nonabelian theory 
such as QCD are expected to be intrinsically different from the case of an ordinary QED 
plasma, because nonabelian gauge fields can directly interact with each other, and because 
charge carriers can have their colors rotated by intervening QCD fields. 

Recently there has been substantial progress on this problem, though it is far from set- 
tled. The growth of QCD plasma instabilities when the soft colored gauge fields are small 
has been well characterized [5, 8]. Arnold and Lenaghan [9] conjectured that exponential 
growth would continue after QCD fields became large. 1+1 dimensional simulations [10, 11] 
confirmed this behavior, but 3+1 dimensional studies [12, 13] showed something different: 
when soft colored gauge fields become large due to instability growth, exponential growth 
shuts off and the energy in soft gauge fields instead grows linearly with time. Recent partial 
attempts to explore the consequences of instabilities for thermalization may be found in 
Refs. [14]. 

The goal of this paper is to further investigate the behavior of an anisotropic quark-gluon 
plasma, after soft gauge fields have become large. We will show that the soft gauge fields 
with wave-numbers which should make them unstable instead enter a dynamical quasi- 
steady state, gaining energy from the instability but losing energy, via their nonabelian 
interactions, to more ultraviolet field excitations. The energy released into infrared gauge 
fields by the instability then cascades towards higher wave-number gauge field modes, with 
occupation number f(k) oc k~ v up to a time- dependent cutoff k max (t). Through simulations, 
we measure the spectral index to be v ~ 2. 

II. SETUP AND FORMALISM 

We will carry out a numerical study of nonabelian, classical Yang-Mills theory coupled to 
an anisotropic bath of high momentum particles. The numerical tools were already presented 
in detail in [12], and this paper is a continuation of that work. To make our presentation 
more self-contained, we will briefly review the setup and formalism here. 

Our treatment is founded on two assumptions about the system under study. The first 
is that the coupling a s <C 1. The second is that there is a separation of scales between the 
momentum of the typical excitation, p, and the screening scale k p. The screening scale 
is set by the momentum p and number density n of the typical excitations as k 2 ~ g 2 n/p. 
The number density is given by n = f p f p ~ p 3 f p , where f p is the occupancy of momentum 

state p and f p is its angular average. Therefore there is a scale separation k <^ p provided 
that the angular-averaged occupancy f p of typical excitations is <C l/g 2 - This is essentially 
a diluteness condition on the typical excitations in the system. 

In the context of heavy ion collisions, the condition a s C 1 is formally valid for the 
collision of extremely large and high-energy heavy ions [6, 15]. The initial angular-averaged 
occupancy is at most f p ~ 1/g 2 , which is the saturation limit [6, 15] 1 (corresponding to 



If the initial occupancies were 3> 1/g 2 , then the magnetic screening length they would provide would be 



3 



number and energy densities n ~ p 3 /g 2 and e ~ p 4 /g 2 )- The density n of the original 
particles, and so f p ~ n/p 3 , will fall with time as the system expands. So the scale separation 
condition f p <C 1/g 2 will be satisfied at least for r ^> t (up until the much later time when 
the particles of momentum p lose their energy and the plasma thermalizes — see, for example, 
the discussion in the original bottom- up scenario of Ref. [6]). 

If we have a clean separation between the screening scale k and typical momenta p, then 
we can describe the physics on length scales x ~ 1/k in terms of (i) classical fields with wave 
vectors ~ k and (ii) a distribution of classical particles (possessing well defined momenta and 
positions) representing excitations with momenta ~ p. This is a Vlasov equation treatment. 
With a little more work [12, 16, 17], one may work in terms of a classical field variable 
A^x), a variable W a (v, x), which represents the net (adjoint) color of all particles moving in 
direction v at point x, and a background colorless particle density with angular distribution 
fl(v, x) and polarizability characterized by a screening mass m 2 , 

with the sum over spin and particle type (including anti-particles). 2 Here t r is a color group 
factor [defined in terms of color generators T by tr(T a T b ) = t r S ab and equal to 3 for gluons]. 
For an isotropic medium, m 2 = \rn{~y = \u) 2 x , where m D is the Debye mass and u p \ is the 
plasma frequency. 

We argued in [5] that the intrinsic length and time scales for the instability to develop 
are short compared to the length and time scales on which the heavy ion system as a whole 
evolves. Therefore we will restrict our attention in this (still somewhat exploratory) study 
to a non-expanding system with spatially uniform VL(v). The evolution equations for 
and W are 

D v F^{x) = [v»W(v,x), (2.2) 

J V 

(D t + v- D X )W = m 2 [E- (2v - V w ) + B ■ (v x V„)] Q{v) . (2.3) 

The first equation is the Yang-Mills field equation, with W(v) giving rise to a current. The 
second equation, derived in Ref. [12], shows how electric and magnetic fields can polarize 
the colorless distribution of particles to create a net color moving in each direction. The 
dynamics of the soft fields is equivalent to that of hard-loop effective theory [18]. 

For a fully nonperturbative study of this system, we implement these equations on a 
lattice. The treatment of the classical, Minkowski space gauge fields on the lattice is the 
standard one [19]. The representation of the W fields is developed in [12, 17] and consists 
of replacing first derivatives with forward-backwards (covariantized) differences in Eq. (2.3), 
and expanding the space of velocities v in spherical harmonics, truncated at a finite but 
large £ max . In this work we make one additional modification which improves approach to 
the £ max — > oo limit and so reduces the memory and time requirements for simulations: We 



larger than the imputed typical momentum. However, the concept of momentum for an excitation of a 
gauge dependent field is only well defined and robust under gauge fixing prescriptions to an accuracy set 
by the magnetic screening length. These ideas are central to the idea of the saturation scale. However 
this argument does not depend on the saturation scenario as such. 
2 In previous work we have referred to m as moo, but here we will drop the subscript. 
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(2.1) 



apply a weak damping term to the large I spherical harmonic components of the W field. 
The reason this is an improvement is that very large i excitations of the W field tend to 
cascade to still larger £, due to the v ■ D term in Eq. (2. 3). 3 With a finite £ max cutoff, some 
of this excitation energy "bounces off" the £ max cutoff and re-enters the infrared. This effect 
becomes smaller as £ max is raised, but can be largely eliminated with the damping term we 
add. This point is discussed and justified in greater detail in appendix A. 

The goal of this paper is to understand the behavior of a would-be unstable system after 
the unstable modes have already grown to have nonperturbative occupancy. In the context of 
a heavy ion collision, this is probably the only relevant question; since the exponential growth 
rate is larger than the inverse system age for all times after the formation time of the plasma, 
there has always been time for unstable gauge field modes to grow to nonperturbative size. 
Therefore we will begin our system with initial conditions for which the infrared (IR) fields 
are already nonperturbatively large. Again, since this is a qualitative exploratory work, 
we will work in SU(2) rather than SU(3) gauge theory and we will only consider a single 
anisotropic particle distribution, the same as the one considered in Ref. [12]. 

In addition to our focus on large rather than small initial conditions, there are slight 
technical differences in our choice of initial conditions compared to those of Ref. [12]. We 
start the system with vanishing E and W fields and with the A field selected from the ther- 
mal ensemble at temperature T = 2m/ g 2 , but then field-smeared to a scale l/2m. 4 This 
is a gauge-invariant procedure which corresponds perturbatively to multiplying a thermal 
spectrum for A(k) by exp(— k 2 /Am 2 ). We will explain the field smearing procedure below be- 
cause it also constitutes one of our best measurables for understanding, in a gauge-invariant 
fashion, what corresponds to infrared phenomena and what to ultraviolet. 

III. RESULTS: BEHAVIOR OF THE NONABELIAN CASCADE 

Our first conclusion, already presented in our previous paper [12], is that the energy 
in soft electromagnetic fields grows linearly with time. This is displayed in Fig. 2, which 
shows that the linear behavior is common to electric and magnetic fields and is robust to 
changes in the lattice spacing. However, this result leaves it unclear whether or not this 
represents continued growth of the soft unstable fields, whether these fields "abelianize" 
[9], and whether they retain long time scale coherence. We now present evidence that the 
answer to all three questions is, "no." 



3 The v ■ D term mixes neighboring £'s because (£m\v\£'m'} has non-zero cases when £' =£ ± 1. 

4 In three dimensions, there should be no important qualitative difference between using a smeared thermal 
distribution for the initial conditions, as here, or a smeared Gaussian noise distribution, as in Ref. [12]. 
We have checked this for a variety of properties, such as the rate of linear growth of energy with time. The 
reason we use a different procedure here than in Ref. [12] is inessential, having to do with the development 
of our code. 
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FIG. 2: (color online) Magnetic energy density J d?xB 2 /2V and electric energy density / d 3 xE 2 /2V 
as a function of time, for three lattice spacings of am = 0.2, 0.288, and 0.4. After a brief transient 
owing to our choice of initial conditions, the electromagnetic energy rises linearly with time. 

A. Chern-Simons diffusion and IR dynamics 

First, consider the behavior of Chern-Simons number, 

N cs (t) - N cs (0) = J dt 1 J d 3 x£^E a (x,t') ■ B a (x,t') . (3.1) 

Chern-Simons number is useful because it characterizes nonperturbative physics. In an 
abelian theory, or a nonabelian theory where the fields are weak, Acs may fluctuate about 
zero but cannot drift away from zero permanently. That is because permanent change re- 
quires topology change (the Minkowski version of instantons). Therefore, ignoring small 
fluctuations, the time evolution of Acs is purely indicative of the dynamics of nonperturba- 
tively large fields. Fully topological algorithms for tracking Acs i n a real-time gauge field 
evolution already exist. We use the one from Ref. [20], which is a modification of a technique 
due to Ambj0rn and Krasnitz [21]. 

A sample Chern-Simons number trajectory, from the am = 0.2 evolution shown in Fig. 
2, is shown on the left in Fig. 3. Chern-Simons number is changing by large amounts 
(indicating nonperturbatively large fields), and in a chaotic fashion (indicating that the 
fields have no long time-scale stability). This means that the soft fields remain large, but 
evolve dynamically, changing on a time scale set by the scale mT x . To see this better, 
consider the auto-correlator (Ncs(t + At)Acs(^)), shown on the right in that figure. This 
correlator indicates over what time scale coherent changes to the gauge fields occur. The 
figure presented represents an average over time and over 10 independent initial conditions 
drawn from the same ensemble, using a 50 3 lattice with am = .288 and £ max = 15. The 
errors were determined by the jackknife method. 
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FIG. 3: Left: typical iVcs evolution in the linear regime. Right: the autocorrelator iVcs-WcSj or 
f x ((E ■ B)[0,t](E ■ B)[x,t + At]), as a function of lag time At. Chern-Simons number changes 
randomly with coherence time scale of order m" 1 . The measurements were made on spatially 
smeared copies of the fields, as described in [20], with smearing extent r = 0.16/m 2 . (/e/fc72 3 
lattice, am = 0.2, £ mSLX = 15; right:50 3 lattice, am = 0.288, £ mSLX = 15) 

Although we have not shown it in the figure, we have also checked that the evolution 
of iVcs does not speed up or slow down during the course of the linear rise in magnetic 
energy. For instance, comparing the first and second halves of the evolutions used to make 
Fig. 3, the Chern-Simons number diffusion rate (sphaleron rate) is consistent between the 
two halves to within 10% error bars, and the full auto-correlator is also consistent within 
errors. What this tells us is that the IR fields truly are nonperturbative, that they evolve 
with a characteristic time scale of order 1/m, and that the time scale and the size of the 
nonperturbative fields is constant throughout the linearly growing regime. 



B. Coulomb gauge spectra 

If IR fields are not growing, then the linear growth of magnetic energy must reflect 
growth of higher-momentum modes of the soft gauge field. To clarify the situation, it would 
be useful to know the power spectrum of the soft gauge fields as a function of wavenumber 
k and time t. We will start with a direct but gauge-dependent measurement of the power 
spectrum in Coulomb gauge. Later we will discuss how the spectrum can be accessed 
indirectly through gauge-invariant measurement involving smearing (also known as cooling) 
of the field configurations. 

Our picture (to be supported by data below) is that, during the linear growth phase of 
the total energy in soft fields, the soft fields consist of (i) a non-perturbative IR component 
plus (ii) a perturbative component in the form of higher-momentum plasmons. We would 
like to know the distribution function f(k) of these plasmons as a function of k. We fix 
Coulomb gauge using the standard algorithm [22] adopted to the real-time case [23] . Then 
we (i) extract the Fourier spectrum of A and E, (ii) evaluate the two point function, 
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averaging over k vectors in narrow blocks of |fc 2 |, 5 and (iii) define the distribution function 
as determined by A and as determined by E through 6 



fA(k) 



N doi V 



f E (k) 



N doi kV 



(E 2 (k)). 



(3.2) 



Here, V is the total spatial volume and 



iVdof = 6 



(3.3) 



accounts for the two transverse polarization states and the 3 = N% — 1 adjoint color states 
in SU(2) gauge theory. 

On scales where the fields are nonperturbatively large, these distribution functions are 
difficult to interpret, depending somewhat on the gauge fixing procedure. We do not expect 
/a to equal indeed, if the dominant fields are slowly evolving or unstable magnetic fields, 
we expect ] a > I'e, perhaps by a large margin. On scales, presumably at larger k, where the 
fields are perturbative, we should see Ja — Je if the relevant degrees of freedom are behaving 
as plasmons with k > m. Therefore, the ratio /a//e serves as a diagnostic of whether the 
physics is nonperturbative and whether the degrees of freedom are primarily independently 
propagating plasmons or something else, such as the magnetic fields associated with hard 
particle currents. 

To test this, we evolved the system for a time of mt = 400 in a 64 3 box with £ max = 15 and 
lattice spacing am = 0.25, corresponding to a (large) physical volume of (16/m) 3 , tracking 
the distribution functions after the initial transient had died and the magnetic field energy 
was undergoing linear growth. The occupancy after initial transients is displayed on the 
left in Fig. 4, showing that, as expected, the IR has nonperturbatively large fields, while 
fields are perturbative at larger wave number. On the right in the figure, we show the 
time development. Each curve is time-averaged over an interval of At = 12.5/m, and the 
central times of consecutive curves are spaced apart by 25/m. The IR occupancies remain 
nonperturbative but with stable amplitude, while the UV occupancy increases. At any k, 
the occupancy rises and eventually saturates; the saturation point progresses to larger k. 
This looks like a momentum-space (kinetic) cascade. 

Cascades of energy from the infrared to the ultraviolet are familiar from turbulence in 
hydrodynamics and from many other physical systems. There are also weakly coupled 
examples such as weak plasma turbulence in traditional plasma physics [24], and theoretical 
studies of post-inflationary thermalization in the early universe [25]. Such cascades typically 
lead to a steady-state, power law distribution for the power spectrum f(k), usually referred 
to as a Kolmogorov spectrum (in honor of the application to hydrodynamic turbulence). 
Different microphysics leads to different powers of k. A thermal spectrum is f{k) oc k~ l . 
Cascades in scalar field theories during "preheating" after inflation, for example, typically 



5 technically, we use the lattice fc 2 , ^^4/a 2 ) sin 2 (A: i a/2). 

6 Consider the total energy of a gas of weakly-interacting, high-momentum (k 3> m) plasmons. This could 
be written in terms of / as iVdof^ J k ^k.fk — J k kfk or in terms of the fields as NdofV J k \{E k + B k ). 
Since E k ~ B\ for such plasmons, and B k ~ k 2 A\ in Coulomb gauge, we can also write the energy as 
j k E k or J k k 2 A k . Comparison of these expressions leads to the identification (3.2). 
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FIG. 4: (color online) Soft gauge field power spectra: left, initial; right, as a function of time. The 
IR fields are in a quasi-steady state, and the energy cascades towards more UV modes. (64 3 latt ice, 
am = 0.25, £ max = 15) 

display a power spectrum with various power laws at different stages, such as / oc k~ 3 ^ 2 
and fc -5 / 3 [25]. Obviously we do not expect power behavior at values of k where the field is 
nonperturbative; indeed, it is not even clear whether to use /a or in this region. However, 
we do expect power behavior for k large enough that /a — /e, but small enough and at late 
enough times that f(k) has become nearly time independent. 

The first figure in Fig. 5 repeats the righthand figure from Fig. 4 and shows that the 
cascade region is well fitted by a power law with / oc k~ 2 . Unfortunately, the most ultraviolet 
wave numbers involved in the cascade are already at large enough k that lattice spacing 
effects may be a concern. 7 As a check on the robustness of the result, we performed a 
second evolution, also in a 64 3 box, with the same choice of Q(v), but with am = 1/6 
rather than 1/4, and going out to a time of mt = 800. The smaller spacing means that the 
physical volume was somewhat smaller. Nevertheless, it was large enough: we have checked 
agreement within errors of the Nqs diffusion coefficients, and close agreement in the energy 
growth rates. In any case the physics of the cascade is presumably more ultraviolet than 
the physics of the instability, and should not show severe volume sensitivity. The power 
spectrum from this evolution is shown on the right in Fig. 5, and is in very good agreement 
with the larger volume figure. The line superposed is drawn to guide the eye and is not 
actually a fit; it is precisely the same line in each figure. 

We have found power-law fall-off f(k) = k~ v with spectral index v ~ 2. To get a crude 
idea of the error in our determination of v, we show a more detailed view in Fig. 6 of the 
late-time distribution on the finer lattice. We show least-square fits of various power laws 
[v = 1.6, 1.8, 2.0, 2.2, 2.4) to the data points 8 in the range 3 < k/m < 8. The v = 2 fit has 
the longest span of agreement, and we take our final result to be v = 2 ± 0.2. 



7 For example, there is a small ripple in the first figure of Fig. 4 at k/m = 8. For am — 0.25, this is the 
lowest k 2 value where the lattice group velocity can vanish, corresponding to a Van Hove singularity. 

8 We fit the /e data. There is only a slight difference between /e and /a, which is at the IR end of the 
range chosen. The reason for slightly preferring /e over /a will become obvious in Sec. Ill C below. 
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FIG. 5: (color online) Left: power spectrum as in Fig. 4, with a slope k~ 2 power law superposed to 
guide the eye. Right: same, but using a finer lattice spacing am = 1/6 and twice as much physical 
time. 




FIG. 6: (color online) The power-law region of the late-time Coulomb gauge distribution / shown 
fitted with lines corresponding to powers v = 1.6, 1.8, 2.0, 2.2, and 2.4. The distribution functions 
are from the long-time, finer-lattice simulation shown in the right-hand plot of Fig. 5, averaged 
over 700 < rat < 800. Each data point represents the center of one of the bins in k 2 used to 
construct averages in (3.2). (64 3 lattice, am = 1/6, ^ max = 15) 

C. Gauge-invariant cooling 

We argued that the Coulomb gauge spectrum of Fig. 5 should be trustworthy away from 
the IR because the higher-momentum components of the field are perturbative (as can 
be seen from the figures by the drop of occupation number with increasing momentum). 
However, in order to be sure that results are not artifacts of gauge fixing, it is usually 
preferable to investigate gauge theories with gauge-invariant observables. It is possible 
to probe aspects of power spectra in a gauge-invariant way by calculating the energy of 
smeared (also known as cooled) gauge fields. Smearing is a gauge-covariant process which 
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is a function of a parameter r known as the smearing depth. Define A(t, x; r=0) to be 
the actual gauge-field configuration A(t, x) at a given physical time t. Then evolve in r 
according to 

dA_ d(B 2 /2) 

l^-—dA—- DjF ^ (3 ' 4) 
Here D is the covariant derivative and F is the field strength, also using the smeared field 
A(t,x;r). This is a gauge-invariant procedure which has a straightforward lattice imple- 
mentation. Such smearings have a long history in lattice gauge theory studies; for instance 
they are also used extensively in our technique for measuring N cs . Other fields can also be 
smeared. For instance, we define a smeared W field through, 

d T W(x, v; t) = D 2 W(x, v, r) , (3.5) 

where again D 2 is the covariant derivative using the smeared gauge field at the same smearing 
depth r. Note that we do not introduce smearing into the dynamical evolution of the fields 
in time (£); we only use smearing for making measurements. To answer the question "what 
is the smeared B 2 at time t?" , we make a copy of the fields at time t, apply smearing, and 
measure B 2 on it. 

Smeared fields are good at telling us whether the energy going into soft electromagnetic 
fields is appearing in very long wavelength, nonperturbative fields, or in plasmons with 
larger wave number. To study this, we consider the magnetic field energy density B 2 /2 
function of time and of smearing depth, shown in Fig. 7. Very roughly, smearing to a depth 
r eliminates fields with k 2 > (2r) _1 and leaves fields with smaller k. But near k 2 ~ (2r) -1 , 
the fields are only partially removed, so the real story is slightly more complicated; smearing 
is similar to a Laplace transform of the power spectrum with k 2 playing the role of time. The 
figure shows that the infrared energy is stable through the evolution. (It bounces around, 
which shows that the fields are evolving and that it is a small number of degrees of freedom 
contributing to the infrared energy. This bouncing would average out in larger volumes 
because of incoherent averaging.) The total energy (r=0) rises linearly. 

Based on the cascade picture of Fig. 5, what behavior should we have expected for the 
intermediate case of moderate r? For any fixed smearing r, the energy should eventually stop 
growing once the modes with k 2 < r have grown to reach their steady-state distribution 
in the cascade. The smaller r, the longer it should take to reach the steady state. The 
m 2 T = 4/256 curve in Fig. 7 is a good example of an initial rise in energy that then tapers 
off and is plausibly approaching a steady-state value. It is unclear from our data whether 
the smaller (non-zero) r curves will eventually reach steady-state values. The problem is 
that cooling does not select out a single k but gives a superposition in the form of a Laplace 
transform. Note, for instance, that the k~ 2 behavior of the Coulomb spectrum of Fig. 4 is 
limited to the relatively narrow range of 3 < k/m < 5. With infinite computing resources, 
one could run long enough, on fine enough lattices, to extend this region over a huge range 
of k, and then one would expect to see the predicted behavior in Fig. 7 for a wide range of 

T. 

But there is a way to check that our Coulomb gauge results are trustworthy. We 
can Laplace transform the Coulomb gauge spectrum and see if it agrees with the (gauge- 
invariant) smeared measurements of Fig. 7. Perturbatively, the smeared magnetic energy 
density should be related to the distributions f(k) by 

W = ^fJ ^kf(k)e- 2 ^. (3.6) 
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FIG. 7: Magnetic energy density as a function of time and smearing depth. (64 3 lattice, am = 0.25, 

^max — 15) 

Note that the growth of d 3 k x k ~ k 3 dk with k typically more than compensates for the 
fall of f(k) with k in Fig. 5, so that these integrals are dominated by k 2 ~ r _1 at late times 
t. We replot Fig. 7 in Fig. 8, superposed with dashed lines that show the results of (3.6) 
with / = /e. The results are very close until one cools deep into the infrared. (We also 
show one sample curve based on / A , which gives slightly less accurate results than f E . The 
/a have a more significant IR contribution than /e, as can be seen in Fig. 4.) We conclude 
that Coulomb-gauge distributions provide a reasonably accurate description of the physics 
of the cascade far from the IR, supporting our conclusions based on Figs. 4 and 5. 

The instability preferentially excites gauge fields with k vector along the z axis (the 
axis about which the particle momentum distribution is oblate) [5, 8]. Such modes have 
primarily transverse magnetic fields, so B 2 <C B 2 +B 2 for the fields excited by the instability. 
Therefore we might expect this behavior of the infrared gauge fields. If the higher momentum 
fields represent a nearly thermalized bath, they will be close to isotropic and the ratio 
B 2 z /(Bl+By) will be 1/2. But if they scatter predominantly off the IR fields, they may also 
carry a momentum space anisotropy. To study this, Fig. 9 shows the ratio B 2 /(B 2 +B 2 ) as 
a function of smearing, for the same lattice parameters as in Fig. 7. Indeed, the soft fields 
are anisotropic as expected. The unsmeared fields are (on average) less so, but still have a 
definite anisotropy along the z axis. 

The picture that has emerged is that there are soft, anisotropic, nonperturbatively large 
gauge fields with higher wave-number plasmons superposed. The size of the soft gauge fields 
fluctuates about a steady mean, and the plasmons become more numerous, populating higher 
and higher wave numbers. One way to check whether the interpretation of the high wave 
number excitations as plasmons is correct, is to look at the W fields. We will concentrate 
on the £ = 1 component of W, which is the same as the particle current up to a factor: 
j 2 = (1/3) ^2 m W 2 m [12]. A plasmon with k ^> m carries almost all its energy in E 2 and 
B 2 , roughly equipartitioned, and only a subdominant amount in currents. Therefore, if 
the energy growth really represents a growing number of plasmons with k > m, then (j 2 ) 
should grow slowly if at all, and should be IR dominated. Fig. 10 shows our measurements, 
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FIG. 8: (color online) The solid lines are the same as Fig. 7. The dashed lines represent the 
corresponding results extracted from the Coulomb-gauge spectra /e of Fig. 4. The single dotted 
line at the top shows a similar extraction of the unsmeared (t=0) curve from /a nstead of /e- The 
noisy difference with the corresponding /e curve is an IR effect, and this difference remains until 
one cools substantially (not shown). (64 3 lattice, am = 0.25, £ max = 15) 
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FIG. 9: The ratio B^/(B^.+By) as a function of time and smearing depth. (64 3 lattice, am = 0.25, 

^max — 15) 



and the lack of growth with time is clear. 9 [Following Ref. [12], we have normalized the 



9 In contrast, a slight growth of current can be seen at very late times in Fig. 11 of Ref. [12]. This growth 
appears to be a late-time artifact of the undamped treatment of ^ ma x in that reference, as we discuss in 
the appendix. 
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FIG. 10: Current squared as a function of time and cooling depth. The current remains predomi- 
nantly infrared, and shows large fluctuations about a nearly flat trend through the evolution. (64 3 
lattice, (ITfl — 0.25, traax — 15) 
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FIG. 11: The ratio of j 2 to j 2 + jy for the various cooling depths of Fig. 10. Note the scale of the 
vertical axis. (64 3 lattice, am = 0.25, £ max = 15) 

curves by plotting (3/4m 2 )j 2 = (l/4m 2 ) W 2 m , which is the contribution of the £=1 
components of W to what would be a conserved energy (E 2 + B 2 )/2 + (l/4m 2 ) ^2e m W 2 m 
for isotropic systems.] j 2 falls more slowly with cooling depth than the magnetic energy 
of Fig. 7, indicating less power in the ultra-violet. We also find, in Fig. 11, that j 2 ~ 
0.13(j 2 +j 2 ), nearly independent of time and smearing depth. A small ratio of j z to j x 
and j y is expected if the currents are primarily associated with long wavelength, transverse 
modes in the directions which are perturbatively unstable. 
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IV. DISCUSSION AND CONCLUSIONS 



Putting our numerical results together, the physical picture which emerges is the follow- 
ing. Plasma instabilities drive IR gauge fields with k z ^> k± to grow. Nonperturbatively 
strong interactions between these soft field modes remove energy from these unstable modes, 
moving it instead into less-IR gauge field modes. The size of the soft nonperturbative fields 
reaches a quasi-steady state; if it grows larger, the nonperturbative physics removing en- 
ergy gets more efficient, and if it gets smaller, the instability drives it back up. The energy 
absorbed via the instability from the hard particles thereby powers a cascade of soft gauge 
field excitation energy towards the ultraviolet. The cascade has power spectrum f(k) oc k~ v 
with v ~ 2. This is not a thermal spectrum, as also evidenced by the failure of B z /(B 2 +By) 
to approach \. We give a theoretical explanation of the value v = 2 in another work [26]. 

These same instability-powered energy cascades should appear in the early stages of 
arbitrarily high-energy heavy ion collisions. These cascades will only be a temporary feature 
along the path to thermalization and will disappear by the time the plasma is finally fully 
thermalized. We leave to future work the complete integration of the physics of instabilities 
into quark-gluon plasma thermalization at arbitrarily high energies. 
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APPENDIX A: IMPROVING SIMULATIONS BY DAMPING HIGH I MODES 

In this appendix we argue that the large £ max limit is achieved more quickly by applying 
weak damping on large £ modes than by not doing so, and we present numerical evidence 
that this is the case. 

At the perturbative level, the behavior of the isotropic version of the A and W field 
system has been investigated by Bodeker, Moore, and Rummukainen [17]. The correct 
analytic structure of the gauge field propagator in the presence of hard thermal loops (the 
infinite £ max , isotropic theory) is that there should be a propagating "plasmon" pole at a 
frequency > k, and a cut in the spectral weight for all frequencies |a>| < k. Physically, 
the cut reflects Landau damping. It means that most of the energy of a long wavelength 
magnetic field should be absorbed by the particle degrees of freedom represented by the 
W fields, never to return. However, the finite £ max system is a non-dissipative Hamiltonian 
system. Therefore it is not possible for it to contain cuts in the (leading order) propagator. 
Instead, the region which should contain the cut contains a series of poles and zeros in the 
propagator, illustrated in Fig. 12, which is taken from Ref. [17]. 

What this means is that the excitation energy present in the magnetic field, which is 
supposed to be Landau damped away, instead appears only in periodic oscillatory modes. 
Most of the energy is stored in the W fields, as is supposed to happen under Landau damping. 
However, a fraction of it periodically reappears as magnetic field energy, an effect which is 
unphysical. The size of this effect is determined by the number of W field modes, (£ max + 1) 2 , 
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HTL A 1 k = OAm D l m = 10 Poles of propagator / max = 10 




FIG. 12: Left: The inverse propagator with £ max = 10, plotted against u/ijid with fixed k = 0.4 mo- 
Right: The positive frequency poles of the propagator at £ max = 10. In these figures, one can 
clearly see the development of the cut in the interval —k<L0<k, and the two plasmon poles at 
uj 2 « m 2 D /3 + 6k 2 /5. 



which sets the heat capacity of these modes to absorb and store excitation energy. The larger 
^max is, the more effectively the W fields can retain the energy rather than feeding it back 
to the magnetic fields. 

The problem in the current context is that, in the anisotropic system, quite large amounts 
of W field excitation are generated. The derivative term v ■ D W in Eq. (2.3) should cause 
this excitation energy to migrate to higher and higher £ values. This is shown in the left- 
hand plot of Fig. 13, which shows how J2 m W 2 m varies with time and with £, if we start 
with large magnetic fields but with no initial excitation in the W fields. However, with a 
finite £ max cutoff, the excitation cannot migrate to arbitrarily large £, as it should; some of 
it instead moves back into lower £ and re-enters the gauge fields and small £ W fields. This 
causes a fake increase in the energy of these fields, which becomes more severe as £ max is 
decreased. 

Our approach to solve this problem is to assume that any excitation energy which reaches 
very large £ values should continue to arbitrarily large £ and be lost to the low £ system. 
We can make this happen artificially by applying a damping term to the highest £ modes, 
modifying the W equation of motion via 

= (previous) - -yW em O(£ - £ damp ) . (Al) 

That is, we add a term which causes exponential shrinkage with rate 7 in all £ at or above 
a cutoff fdamp- The damping should only affect the high £ modes, so we choose £dam P = 
(2/3)£ max . Note that the damping we apply is gauge invariant: so long as £dam P > 1, the 
current remains exactly conserved, because the current is determined by the £ = 1 modes 
and the charge density by Wqq- We also have to choose a value for 7. If the value is too 
large, then the W fields with £ > £d amp are effectively forced to be zero, which is equivalent 
to lowering £ max to £da,m P - Therefore we should choose 7 < m, since m is the intrinsic scale 
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FIG. 13: Time and £ dependence of W^ m the amount of excitation in the W field. Each figure 
plots J2 m against £ at a series of time snapshots, for £ max = 36. The left figure is without 
damping; the right figure has damping on £ > 24 modes. In each case, the time between successive 
lines is 0.57/m. The two spikes at small £ are the angular scales driven by the instability. The 
excitation introduced here cascades to larger £, eventually bouncing off the £ max limit in the left 
figure but being absorbed at large £ by the damping in the right figure. (50 3 lattice, am = 0.288) 



of the dynamics. We should also make sure we pick 7 > m/£ niax ; otherwise the damping is 
too slow to do anything, as excitation energy can get from £dam P to £ max , reflect, and go back 
below fdamp before being damped away. With this in mind, we have chosen 7 = m/\/£ max . 

Of course, it remains to test whether this procedure of damping large I modes really 
gives the same behavior as choosing an extremely large value for £ max (which is numerically 
impractical, since the number of degrees of freedom rises as (l+^ ma x) 2 )- To check, we have 
performed evolutions with the same initial random seed and other parameters but with 
different values of £ max , and either with or without the damping term added. The results 
are shown in Fig. 14 for am = .288, V = (14.4/m) 3 , and the same fl(v) used in the main 
body of the paper [12]. On the left, we see magnetic energy growth vs. time. With and 
without damping, the growth rates converge (from opposite sides) towards the same large 
£ max limit. But the damped results converge much faster: with damping, the magnetic 
energy growth does not change between £ max = 15 and £ max = 24, and appears to coincide 
with the very large £ max , undamped behavior. Now return to the flat result of Fig. 10 for 
the time dependence of j 2 . On the right of Fig. 14, we see that this flat behavior is obtained 
only very gradually in the large £ max limit unless damping is implemented, in which case 
it occurs immediately. For investigating current growth, damping is an essential numerical 
tool for practical simulations. 

Naturally, if £ max is too small, we will see errant behavior whether or not we implement 
damping. For the choice of Q(v) used here (i.e. for the degree of hard particle anisotropy 
in our simulations), Fig. 14 shows large deviation of the magnetic energy growth for £ max 
below about 10. We have therefore conservatively used £ max = 15 for the bulk of the studies 
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FIG. 14: (color online) Left: magnetic energy against time at different values of l mSuX , with and 
without W damping. Without damping, the growth rate in B 2 approaches a large ^ max value from 
above, approximately as l/(4nax+l) 2 - With damping, it approaches from below and obtains the 
large £ max value much faster. Right: the same for (j 2 ). The fact that this quantity does not grow 
is obtained immediately with damping, but only approached very slowly in the large ^ max limit 
without damping. (50 3 lattice, am = 0.288) 
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